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Singular perturbations and Lindblad-Kossakowski 

differential equations 

Mazyar Mirrahimi and Pierre Rouchon 



Abstract — We consider an ensemble of quantum systems whose 
average evolution is described by a density matrix, solution of 
a Lindblad-Kossakowski differential equation. We focus on the 
special case where the decoherence is only due to a highly 
unstable excited state and where the spontaneously emitted 
photons are measured by a photo-detector. We propose a sys- 
tematic method to eliminate the fast and asymptotically stable 
dynamics associated to the excited state in order to obtain 
another differential equation for the slow part. We show that 
this slow differential equation is still of Lindblad-Kossakowski 
type, that the decoherence terms and the measured output 
depend explicitly on the amplitudes of quasi-resonant applied 
field, i.e., the control. Beside a rigorous proof of the slow/fast 
(adiabatic) reduction based on singular perturbation theory, 
we also provide a physical interpretation of the result in the 
context of coherence population trapping via dark states and 
decoherence-free subspaces. Numerical simulations illustrate the 
accuracy of the proposed approximation for a 5-level systems. 

Index Terms — Quantum systems, Lindblad-Kossakowski mas- 
ter equation, singular perturbations, optical pumping, coherent 
population trapping, adiabatic approximation. 



L Introduction 

Under the usual assumptions of optical pumping and/or co- 
herent population trapping, the Lindblach-Kossakowski master 
equation describing the dynamics of the density operator 
admits multiple time- scales. In this paper, we are studying 
the fast/slow structure resulting from a separation between 

• the life-times of the excited and unstable states assumed 
to be short. 

• the oscillation periods, assumed to be long, associated to 
the energies of the other stable states and to the Rabi 
pulsation generated by the control, coupling the unstable 
and stable states in a quasi-resonant way. 

Usually, the elimination of the fast dynamics is performed in 
terms of coherence vector gathering in a single column the 
coefficients of the density matrix. In this form, the system 
is not written in a standard form |7, chapter 9] (also called 
Tikhonov normal form) with a clear splitting of the coherence 
vector into two sub- vectors: a fast sub- vector and a slow one. 
Some tedious linear algebra and changes of variables are then 
needed to put the system into the standard form in order to 
perform the adiabatic (quasi-static approximation) elimination 
of the fast dynamics. Moreover with such coherence vector 
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we loose the physical interpretation of the master equation 
in terms of Hamiltonian and jump operators, explained in (61 
chapter 4]. 

The main contribution of this note is to propose a more 
intrinsic elimination of the fast part of the dynamics by 
using only matrix manipulations for systems with a structure 
sketched on figure [T] The main theoretical guide is the 
geometric theory of singularly perturbed differential systems 
initiated in [5] and center manifold techniques to approximate 
the invariant slow manifold Q, O. Such theoretical guides 
have been already used in f?! in the context of reduction of 
kinetics combustion schemes. These guides avoid here the use 
of the coherence vector and provide a slow dynamics that is 
also a Markovian master equation of Lindbald-Kossakowski 
type with a slow Hamiltonian and slow jump operators. This 
slow master equation describes the dynamics of the density 
matrix of the open quantum system that lives in the Hilbert 
space spanned by the stable states. As far as we know, such 
formulation of the slow dynamics is new, even in the physicist 
community, and could be of some interest for the control. In 
particular, the controls appear explicitly in the decoherence 
terms and the output map. 

The note is organized as follows. Section |ll|is devoted to the 
modeling of systems depicted on figure [T] to the three time- 
scales structure and to the elimination (by averaging, i.e., by 
the rotating wave approximation usually used by physicists) 
of the fastest time- scale attached to the transition frequencies 
between the stable and unstable states. The resulting model, 
equation ^ with complex value controls Vik and measured 
output still admits two time- scales, an asymptotically stable 
fast part and a slow part. Extraction of the slow part is the 
object of section [In] where the approximation Theorem [T] is 
proved. For readers not interested by these technical devel- 
opments, we have summarized at the end of this sectio n [lll| 
the main formula for deriving the slow master equation ( |23| ) 
from the original slow/fast one. In section IV we compute 



the slow approximation when the Hamiltonian H corresponds 
to ([5]) and provides physical interpretations in terms of slow 
Hamiltonian and slow jump operator depending directly on the 
control amplitudes I]/.. In section |Vj we compare, numerically 
on a five-level system, the slow/fast master equation with the 
slow one. 

A preliminary version of these results can be found in O. 
The authors thank Guilhem Dubois from LKB for several 
discussions on the physics underlying coherence population 
trapping. 
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II. The three time-scale master equation 

Such master equations typically describe coherent pop- 
ulation trapping when a laser irradiates an (TV + 1) -level 
system |2|, |1|. The system is composed of N (fine or 
hyperfine) ground states {\gk))k=i having energy separations 
in the radio-frequency or microwave region, and an excited 
state |e) coupled to the lower ones by optical transitions at 
frequencies {oOk)^^i (see Figure [l]). Naturally, the decay times 
for the optical coherences are assumed to be much shorter 
than those corresponding to the ground state transitions (here 
assumed to be metastable). 




5fp 



Fig. 1. relevant energy levels, transitions and decoherence rates for the 
considered model ^ 



The quantum Markovian master equation of Lindblad- 
Kossakowski type, modeling the evolution of a statistical 
ensemble of identical systems given by Figure [T] reads (see |6l 
chapter 4] for a tutorial and up-to-date exposure on such 
master equations): 



+ 



y 



^ ^ Tk{2QkpQi - QiQkP - pQlQk) (1) 

k=l 



k=l 



(2) 



where G M is the controlled input (laser field) and y > 
is the measured output (number of photons per time unit 
spontaneously emitted from the excited state |e)). Here the 
Hermitian operators Hq and Hi are, respectively, the free 
Hamiltonian and the interaction Hamiltonian with a coherent 
source of photons u{t) G M: 



N 



- = A|e) (e| + ^A/e \gk) {gk 



H 



N 



the pk being coupling and constant parameters. Moreover, the 
quantum jump operators Qk corresponding to the spontaneous 
emission from the state |e) towards \gk) are given as follows 

Qk = \gk) (e| . 
One easily has the following relations 

QkQi = 0, QlQk = P=\e){e\ Vfc ^ / G {1, • • • , N}. 

(3) 

The transition frequencies, cok = A — A/e (where A and Xk are 
the eigenvalues of ^ corresponding to the energy levels |e) 
and \gk), respectively), are supposed to be much larger than 
the decoherence rates Tk. The control field u{t) is assumed 
to be in the quasi-resonant regime with respect to the natural 
frequencies of the system: 



u{t) 



N 
k=l 



(4) 



where the complex amplitudes ii/e G C are varying slowly and 
where Sk are the small de-tuning frequencies. We have thus 
three different time scales: 

1) The very fast time-scale associated to the optical fre- 
quencies COk- 

2) The fast time- scale associated to the life times of the 
excited state |e), Tk- 

3) The slow time-scale associated to the laser amplitude 
\pkUk\ and to the other atomic transition frequencies 

^ki = ujk - k ^ I. 

We are interested here by the slow time-scale of system ([T]) 
where the control u{t) is given by ^ with the following time- 
scales separation: 



and 



d 



< ^k'\Uk\ 



with k,l,k',k" G {1,...,7V}, A:7^/. 

Elimination of the fastest time-scales is standard. It corre- 
sponds to the so-called rotating wave approximation and can 
be justified by averaging techniques. This is not the object 
of this note and thus we just recall here the resulting master 
equation: 

d 1 ^ 

j^P = -\[H^ + 2 ^ ^ki^QkpQl - QlQkP - pQlQk)^ 

k=l 

(5) 

Calculating the secular terms of ue * ^ ^i^ie* ^ % the effec- 
tive Hamiltonian H is given as follows: 



H 



N 



\gk) {gk \ + ^/c \gk) (e| + 1^^ |e) {gk\ • (6) 



/c=l 



/c=l 



with Vtk = PkUk' Note that the measured output remains 
unchanged 

y = (E j Tr (Pp) = Tfc j Tr (|e) (e| p) . (7) 
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We are led to the following master equation where ^ only appears in first equation defining ^p/. There- 

fore p/ is associated to the fast part of the dynamics and ps 
represents the slow part. The fast part is asymptotically stable 



d 



^^hlgk) {9k\^^k\9k) {e\^Vtl\e) {gk\ , p 



./c=l 



^ because — -^^ — — iPf~^ ^Pf^) defines a negative definite 

£^(2 {e\p\e) \gk) {gk\ - \e) {e\ p - p \e) {e\ ) super-operator on the space of Hermitian operators: 

k=i ^^f_f ^„ ^ — j|2 I II pii2 



(8) 



Tr(-(p^ + PpfP)pf) = -{WpfW + \\PpfP\n- 



/ N \ Moreover its inverse is given by : 

y= K]rJTr(|e)(e|p) ^ 

Vfe=i / X^X--PXP. (14) 



where the Q^'s are the slowly varying complex amplitudes ^^^^ ^^^^^ ^1^^ ^^^^^^^ approximation (|25l 

(controlled inputs), the 4 's are the laser de-tunings and where described in the Appendix lAl Computing the first order terms, 
the two time-scales separation results from: following approximation for pf with respect to p., : 



\5k\,\^k\^Tk> and 
for k,k' e {1,...,N}. 



~l~^k 

dt 



« Tfe'l^^fel Pf = . I' , {PHp, - p,HP) + O(e^). (15) 

Inserting now the equations ([TT]) into the equation ([13]), we 
have: 

III. Slow/fast reduction , 
We can therefore take ik = 1 fc/e where e is a small positive dt h ^l^^jv / /l^^^jv y 



^{l-P)J2^k[H, QkPfQi] (1 - P) 



parameter. Thus we have a master equation with the following 
structure: 

d ^ T 

where ffc's and H (given by ^ for example) have the same ^^^^^ ^^^^ ^^^^ ^ 
orders of magnitude but where e > is a small parameter. 

Define, with P = \e) (e|, QkPs = PsQI = 0. (16) 

pj, — Pp^ pP — PpP Applying now the first order approximation ([T5j, and after 

^ some simple but tedious computations, we have 

We have y E/c=i j 

1 A- t ^ ps(l- P)HPH(1- P)) 

P = Ps^Pf- Y^r.Q.pfQl (11) ^ ^ V 

ErfcQfci7p,i7Q^ + 0(e2). (17) 



7 1 A N 

k=i 4e 



=1 



and therefore p {pf^ps) is a bijective map. This map is a fi^f T/c) 

sort of "change of variables" decoupling the slow part from 

the fast part of the dynamics. Note that, in the slow part, p„ Here, we have in particular appHed ^ as well as the fact 

we have somehow removed the fast dynamics associated to that Q^P = Qk- Continuing the computations, we have 
the optical state |e). Indeed, this change of variable leads to a ^ ^ 

standard form: —ps = —-[H^ps]-\- 

dt a 



d 



[T.k=l^k 



N 

{pf + PpfP) 2e Tk {2QkPsQi - QkQkPs - PsQlQk) (18) 



h 



{P[H,p] + [H,p]P - P[H,p]P), (12) ^^^^^ ^^^^ defined 



ih^Os = (1 - P)[H, P](l-P) H = {1- P)H{1 - P) (19) 



df 



1 and 

(^N .Y.'^kQk[H,p]Q^. (13) ^ , ^ _. (l-P)QuH(l-P). (20) 
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Note that 



Denoting hy Sps = Ps — Ps, we have 



QlQk 



Al-P)HPH{l-P), 



which, in particular, allows us passing from ([TT]) to ([18]). 

The situation is different for the measured output y. We 
have: 



y{t) = ^ ^Tr(Pp) = ^ '-Tr{Ppf) 



Se^Tk [Tr (QkSpsQiSps) -TtIqIq^Sp, 

+ Tr (O(e^)^p. 

This, together with Cauchy-Schwartz inequality, implies 



-2z 



Tr{P{PHps-psHP))^0{e). 



Tr (Sps\t)) <Tr(sps\o)) ^eL Tr^ 

rt 



^Ps (r) 



But Tr {P{PHps - psHP)) = 0. We should therefore con- 
sider the second order terms otherwise the first order approxi- 
mation yields y = 0. Using the Appendix [A] simple but tedious 
computations end up by the following natural approximation: 

N 

^(t) =4e(^rfc)Tr(Pp,) +0(6^), (21) 



[ Tr^ 

Jo 



dr 



where L and C are positive constants. Note that, Sps being 
definite positive, we have 



Tr2 



k=l 



where we have defined 



Therefore, noting ^ = W Tr 



4 


< Tr 


^ 2 

Sps (r) . 


h 


" ^ 2 

^Ps (r) 


, we have 



P — QkQk 



n^Eti'^k) Jo Jo 



In order to derive ( [2T] ), we only need to apply ([26]) with the 
appropriate values of the functions given in ([12]) and ([TSj and 
the inverse map given in ( [T4j ). 

We can therefore prove the following theorem: 

Theorem 1. Consider p the solution of the Markovian master 
equation ^ and ps the solution of the slow master equa- 
tion ( [TS] ) with ( [T9l) and (|2Qj). Assume for the initi al states 
|p(0)-p,(0)| =V7>-((p(0)-p.(0))(p(0)-p,(0))) = 0(e). 
Then 



Denoting C, = ^(t) + ^e, some simple computations lead to 



C\t) <2C\0) + 2eL [ C^{T)dT 
Jo 



<e{0) + ^e' + '2eLj\\r)dr. 



Applying the Gronwall lemma, we have 

C2 



Cit) < 



r(o) 



2L^ 



„2eLt 



\p{t) - Ps{t)\ = VTrMt)-p,mp{t)-p,m = 0(e) 

on a time scale t ^ 1/e. 

Moreover the output y{t) of the system (given by ^) may 
be written as in ( [2T] ). 

Remark 1. Note that, the approximation of this theorem is 
stronger than the usual one only ensuring an error of order 
0(e) on a finite time T rather than on a time scale oft^ T/e. 
This stronger result is due to the Hamiltonian structure of the 
dominant part of the dynamics. 

Proof: Applying ( [TT] ) and the singular perturbation theory 
of the appendix [a| we have p{t) = ps{t) + 0(e) where 

d ^ I — ^ 

AT 

2e5] Tfc {iQ^PsQl - QlQkPs - PsQlQk) + 0{e'), 



Noting that, by the Theorem's assumption, ^(0) = 0(e), we 
have ({t) = 0(e) on a time scale of t ~ 1/e. As ^(t) = 
C{t) + 0(e), this trivially finishes the proof. ■ 
From a practical point of view, the main result of this section 
is as follows. The correct slow approximation (also called by 
physicists adiabatic approximation) of the system described by 

^ P = -li^^P] 



dt 



E Y {^QkpQk - QkQkP - pQlQk) 



k=l 



with Qk = \gk) (e| and where the F^'s are much larger than 
^ and where the output reads y = '^j.^kT'^ ^Q^QkP^ , is 
given by 

d 



dV 



-\Hs,Ps 



k=l 



2^ -Y~ [^Qs,kpsQlk - Qs,kQs,kps - psQ 



Ps{0) = p{0). 



(22) 



k=l 



kPs Ps^g^k^^^k 

(23) 
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where ps is the density operator associated to the space 
spanned by the |^/c)'s, where the slow Hamiltonian is 

Hs = (l- P)H(l - P) 

and the slow jump operators are 

Qs,k = Qk^{l-P)- 

Here, we have set T = Xl/c^/e and P = \e) {e\. The slow 
approximation of the measured output is still given by the 
standard formula 

n 
k=l 

IV. Physical interpretation 

In this section, we provide a physical interpretation of 
the last section's result for the particular Hamiltonian of the 
system dSl). We get 



Hs = ^Sk \gk) {gk 



and 



\gk){bn\ with \bn)- ^'"''^'^ 



Let us set 



7fe 



the slow master equation reads 



4r, 



d_ 
dt 



Ps 



\gk) {gk\ , Ps 



+ ^ ^(2 {hn\p\hn) \gk) {gk\ - \hn) {hn\ p - p \hn) {hn\ ) 

k 

with ys = {^k^k) {^n\ps\^n)- Thus, whenever all the de- 
tunings 5k vanish, the unitary state is the bright state 
and the vector- space orthogonal to 1 6^7) is a decoherence free 
space since on this sub-space, the Lindbald-Kossakowski terms 
identically vanish and the output y is null. Notice that the 
controls Vtk appear only in the decoherence terms and have 
disappeared form the slow Hamiltonian. 

If we restrict ourselves to the case of a 3-level A-system 
{N = 2), we have Hs = f (|^i) {gi\ - \g2) (^2!), 



1^1) + 



no 



\92) 



is the bright state of the A-system, as in the context of 
the coherent population trapping. As it can be seen easily, 
whenever no de-tuning is admitted (6 = 0), the dark state 



Id) = \b)^ 



\9i) 



I52) 



is the only equilibrium state of the slow dynamics. 



V. Numerical simulations 

Finally let us check the relevance of the adiabatic reduction 
result of the Section [IV] in a simulation. Here, we consider a 
(4+l)-level system given by the following parameters 

^1 = 0.5 52 = 1.2 53 = 0.7 54 = 1.0 

= 1.0 = 1.2 = 1.1 ^4 = 1.3 
Ti = 5.0 r2 = 4.0 Ts = 7.0 = 5.0 (24) 

The simulations of Figure [2] illustrate the output signal derived 
from the reduced slow dynamics, ys, versus the slow/fast 
dynamics, y. The simulation time is taken to be 2.5Ts, where 
Ts, the time scale of the slow system is given by 

Fi + + Fa + r4 



The initial conditions are identical p{0) = Ps{0) = 
T.k=i\9k){9k\ ^ observe, for the slow/fast dynamics, a first 
fast transient corresponding to the relaxation time of the fast 
dynamics for t between and l/F^, i. e. t G [0, 1/4], and then 
two slow transients very similar for both master equations: the 
slow approximation is clearly valid only for time-scale much 
larger than 1/F/e. 



time 

Fig. 2. The output of the reduced slow dynamics ( [23) (solid line) versus the 
output for the (4+l)-level system jsj (dashed line). The parameters are listed 
in {24} . Even if the time-scale separation is not so large with an e around 
^ , such adiabatic approximation captures quite precisely the slow part of the 
dynamics, as stated by theorem [T] 



Conclusion 

We observed that for an ensemble of independent and 
identical quantum systems, and whenever the decoherence 
dynamics due to the measurement is much faster than the other 
dynamics, the adiabatic approximation helps us to find the 
slow dynamics as well as the measurement result with respect 
to the slow dynamics. Note that in this new system, the deco- 
herence term can be removed in a first order approximation. 
We obtain therefore a system of the form 

:Ps = -TY^s.Psh 



dr 



where the control appears linearly in the reduced Hamiltonian 
Hs. This system corresponds to a bilinear system with the 
wavefunctions as state variables. Furthermore, we have access 
to a measurement y given by the reduced slow evolution. We 
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can henceforth consider a control problem with continuous 
measurement associated to this system. 

Notice that, by some simple but tedious computations, one 
can extend the result of this paper to the more general case of 
an + M-level system with TV metastable ground states and 
M highly unstable excited states. 
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Appendix 

This appendix has for goal to remind an approximation 
technique that can be perfectly justified using the geometrical 
tools of singular perturbation and the center manifold Q, O, 

111 

Consider the slow/fast system {x and y are of arbitrary 
dimensions, / and g are regular functions) 

where x and y are respectively the slow and fast states 
(Tikhonov coordinates), all the eigenvalues of the matrix A 
have strictly positive real parts, and e is small strictly positive 
parameter. Therefore the invariant attractive manifold admits 
for equation 

y = eA-^g{x,{))^0{^) (25) 

and the restriction of the dynamics on this slow invariant 
manifold reads 

j^x = f{x,eA-'g{x,Q)) + 0{e^) 

= fix, 0) + e^\^,,o)A-'gix, 0) + 0{e^). 

The Taylor expansion of g can be used to find the higher order 
terms. For example, the second order term in the expansion 
of y is given by: 

y = eA-^g{x,0)^ 

e'A-' (^^\^^^,^A-'g{x,0)-A-'^\^^^^^ 

(26) 

and so on. 



